{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Elasticity\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.geom2d as geom2d\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a rectangle with refinement at corners:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = geom2d.SplineGeometry()\n",
    "pnums = [geo.AddPoint (x,y,maxh=0.01) for x,y in [(0,0), (1,0), (1,0.1), (0,0.1)] ]\n",
    "for p1,p2,bc in [(0,1,\"bot\"), (1,2,\"right\"), (2,3,\"top\"), (3,0,\"left\")]:\n",
    "     geo.Append([\"line\", pnums[p1], pnums[p2]], bc=bc)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cauchy-Green tensor and hyperelastic energy density:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E, nu = 210, 0.2\n",
    "mu  = E / 2 / (1+nu)\n",
    "lam = E * nu / ((1+nu)*(1-2*nu))\n",
    "\n",
    "def C(u):\n",
    "    F = Id(2) + Grad(u)\n",
    "    return F.trans * F\n",
    "\n",
    "def NeoHooke (C):\n",
    "    return 0.5*mu*(Trace(C-Id(2)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "stationary point of total energy:\n",
    "$$\\delta \\int W(C(u)) - f u = 0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor = Parameter(0)\n",
    "force = CoefficientFunction( (0,factor) )\n",
    "\n",
    "fes = H1(mesh, order=4, dirichlet=\"left\", dim=mesh.dim)\n",
    "u  = fes.TrialFunction()\n",
    "\n",
    "a = BilinearForm(fes, symmetric=True)\n",
    "a += Variation(NeoHooke(C(u)).Compile()*dx)\n",
    "a += Variation((-InnerProduct(force,u)).Compile()*dx)\n",
    "\n",
    "u = GridFunction(fes)\n",
    "u.vec[:] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a simple Newton solver, using automatic differentiation for residual and tangential stiffness:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveNewton():\n",
    "    res = u.vec.CreateVector()\n",
    "    \n",
    "    for it in range(10):\n",
    "        print (\"it\", it, \"energy = \", a.Energy(u.vec))\n",
    "        a.Apply(u.vec, res)\n",
    "        a.AssembleLinearization(u.vec)\n",
    "        inv = a.mat.Inverse(fes.FreeDofs() ) \n",
    "        u.vec.data -= inv*res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "factor.Set(0.4)\n",
    "SolveNewton()\n",
    "scene = Draw (C(u)[0,0], mesh, deformation=u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor.Set(factor.Get()+0.4)\n",
    "SolveNewton()\n",
    "scene.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute $2^{nd}$ Piola Kirchhoff stress tensor by symbolic differentiation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_=C(u)\n",
    "sigma = NeoHooke(C_).Diff(C_)\n",
    "\n",
    "Draw (sigma[0,0], mesh, \"Sxx\", deformation=u, min=-10.001, max=10.001)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
